QK Box 12.2

Chadha et al. (2019) set up a completely randomized experiment studying the effects of different soil water holding capacities (WHCs) on various characteristics of the agronomic weed Lactuca serriola. Plants were grown individually in pots allocated to one of four WHCs (100%, 75%, 50% and 25%). There were seven plants for each WHC, although we just use a subset of the data from the 100% treatment group so there was no between-subjects factor. The number of leaves on each of the seven plants in 100% WHC soil was recorded weekly for nine weeks (counts at the start of the experiment, week 0, were omitted). Time was the within-subjects (repeated measures) fixed factor and individual plants were the random subjects. With nine weeks and reasonably linear trends through time for each plant (Figure 12.2), it made sense to treat time as a continuous covariate for analysis. Although week 0 was not included in the analysis, we did not centre time for analysis, so intercepts represent the number of leaves for week 0.

Prickly lettuce or Milk thistle. Mick Keough, CC SA-BY 4.0

**Link to paper **

Chadha, A., Florentine, S., Chauhan, B. S., Long, B. & Jayasundera, M. (2019). Influence of soil moisture regimes on growth, photosynthetic capacity, leaf biochemistry and reproductive capabilities of the invasive agronomic weed; Lactuca serriola. PLoS One, 14, e0218191.

Preliminaries

First, load the required packages (afex, car, lattice, lme4, lmerTest, nlme, VCA, ez, emmeans, Rmisc, MuMIn)

Import chadha data file

chadha <- read.csv("../data/chadha.csv")
chadha

NA

set contrasts from afex

set_sum_contrasts()

setting contr.sum globally: options(contrasts=c('contr.sum', 'contr.poly'))

select 100%WHC

chadha1 <- subset(chadha, treat=="100WHC")

make individual a factor

chadha1$plant <- factor(chadha1$plant)
options(digits=10)

plot slopes - very consistent

xyplot(noleaves~week|plant, type=c("p","r"), auto.key=T, chadha1)

Fit OLS “ancova” model

chadha1.aov <- aov(noleaves~week*plant, data=chadha1)
plot(chadha1.aov)

summary(chadha1.aov)

            Df   Sum Sq  Mean Sq   F value     Pr(>F)    
week         1 4999.050 4999.050 310.55030 < 2.22e-16 ***
plant        6  710.095  118.349   7.35207  1.221e-05 ***
week:plant   6  110.400   18.400   1.14304    0.35199    
Residuals   49  788.772   16.097                         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get Type III SS)

Anova(lm(chadha1.aov), type='III')

Anova Table (Type III tests)

Response: noleaves
               Sum Sq Df   F value  Pr(>F)    
(Intercept) 3567.6042  1 221.62622 < 2e-16 ***
week        4999.0500  1 310.55030 < 2e-16 ***
plant        130.6516  6   1.35272 0.25245    
week:plant   110.4000  6   1.14304 0.35199    
Residuals    788.7722 49                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get correct F-ratio and P value for week (tested against week by plant)

4999.050/18.400

[1] 271.6875

1-pf(4999.050/18.400, 1, 6, lower.tail = TRUE, log.p = FALSE)

[1] 3.178231089e-06

no GG and HF adjustments as week is continuous with 1 df

Get variance components

chadha1.vca <- anovaMM(noleaves~week+(plant)+(week*plant), chadha1)
chadha1.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

      int      week 
16.400794  3.450000 


    [Variance Components]

  Name       DF        SS         MS         VC        %Total    SD       CV[%]    
1 total      15.101789                       26.537662 100       5.151472 15.308619
2 plant      6         710.095238 118.349206 10.401893 39.196719 3.225197 9.584311 
3 week:plant 6         110.4      18.4       0.038377  0.144613  0.1959   0.582156 
4 error      49        788.772222 16.097392  16.097392 60.658668 4.012156 11.922915

Mean: 33.650794 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(chadha1.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

    int    week 
16.4008  3.4500 


    [Variance Components]

  Name       DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
1 total      15.1018                   26.5377 100     5.1515 15.3086        
2 plant      6       710.0952 118.3492 10.4019 39.1967 3.2252 9.5843  77.6065
3 week:plant 6       110.4    18.4     0.0384  0.1446  0.1959 0.5822  0.0343 
4 error      49      788.7722 16.0974  16.0974 60.6587 4.0122 11.9229 10.5766

Mean: 33.6508 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
           Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total       26.5377 14.5063 63.3391       15.9486       54.6613
plant       10.4019 -6.8643 27.6681       -4.0884       24.8922
week:plant   0.0384 -0.3245  0.4013       -0.2662        0.3429
error       16.0974 11.2325 24.9968       11.8901       23.2468

> SD:
-----
           Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total        5.1515  3.8087 7.9586        3.9936        7.3933
plant        3.2252 -2.6200 5.2600       -2.0220        4.9892
week:plant   0.1959 -0.5697 0.6335       -0.5159        0.5856
error        4.0122  3.3515 4.9997        3.4482        4.8215

> CV[%]:
--------
           Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total       15.3086 11.3183 23.6505       11.8677       21.9707
plant        9.5843 -7.7858 15.6313       -6.0087       14.8264
week:plant   0.5822 -1.6929  1.8825       -1.5332        1.7403
error       11.9229  9.9596 14.8575       10.2470       14.3280


95% Confidence Level  
SAS PROC MIXED method used for computing CIs 

Drop interaction

chadha2.aov <- aov(noleaves~week+plant, data=chadha1)
Anova(lm(chadha2.aov), type='III')

Anova Table (Type III tests)

Response: noleaves
               Sum Sq Df   F value     Pr(>F)    
(Intercept) 3567.6042  1 218.22097 < 2.22e-16 ***
week        4999.0500  1 305.77874 < 2.22e-16 ***
plant        710.0952  6   7.23911 9.9979e-06 ***
Residuals    899.1722 55                         
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

get new variance components

chadha2.vca <- anovaMM(noleaves~week+(plant), chadha1)
chadha2.vca

ANOVA-Type Estimation of Mixed Model:
--------------------------------------

    [Fixed Effects]

      int      week 
16.400794  3.450000 


    [Variance Components]

  Name  DF        SS         MS         VC        %Total    SD       CV[%]    
1 total 23.462944                       27.681988 100       5.261368 15.635196
2 plant 6         710.095238 118.349206 11.333402 40.941432 3.366512 10.004257
3 error 55        899.172222 16.348586  16.348586 59.058568 4.043338 12.015581

Mean: 33.650794 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA

VCAinference(chadha2.vca, alpha=0.05, VarVC=TRUE, excludeNeg=FALSE, constrainCI=FALSE)

Inference from Mixed Model Fit
------------------------------

> VCA Result:
-------------

    [Fixed Effects]

    int    week 
16.4008  3.4500 


    [Variance Components]

  Name  DF      SS       MS       VC      %Total  SD     CV[%]   Var(VC)
1 total 23.4629                   27.682  100     5.2614 15.6352        
2 plant 6       710.0952 118.3492 11.3334 40.9414 3.3665 10.0043 57.76  
3 error 55      899.1722 16.3486  16.3486 59.0586 4.0433 12.0156 9.7191 

Mean: 33.6508 (N = 63) 

Experimental Design: balanced  |  Method: ANOVA


> VC:
-----
      Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total  27.6820 16.7947 54.0455       18.1687       48.3233
plant  11.3334 -3.5623 26.2291       -1.1675       23.8343
error  16.3486 11.6201 24.7038       12.2651       23.0805

> SD:
-----
      Estimate  CI LCL CI UCL One-Sided LCL One-Sided UCL
total   5.2614  4.0981 7.3516        4.2625        6.9515
plant   3.3665 -1.8874 5.1214       -1.0805        4.8820
error   4.0433  3.4088 4.9703        3.5022        4.8042

> CV[%]:
--------
      Estimate  CI LCL  CI UCL One-Sided LCL One-Sided UCL
total  15.6352 12.1784 21.8466       12.6668       20.6578
plant  10.0043 -5.6088 15.2194       -3.2109       14.5079
error  12.0156 10.1300 14.7702       10.4073       14.2767


95% Confidence Level  
SAS PROC MIXED method used for computing CIs 

Now mixed effects modelcomparing random slopes and random intercepts using ML

chadha.lmer1 <- lmer(noleaves~week + (week|plant), REML=FALSE, chadha1)

boundary (singular) fit: see help('isSingular')

 # singular fit so set correlation between random effects to zero
chadha.lmer1 <- lmer(noleaves~week + (week||plant), REML=FALSE, chadha1)
chadha.lmer2 <- lmer(noleaves~week + (1|plant), REML=FALSE, chadha1)

Compare models

anova(chadha.lmer1, chadha.lmer2)

Data: chadha1
Models:
chadha.lmer2: noleaves ~ week + (1 | plant)
chadha.lmer1: noleaves ~ week + (week || plant)
             npar       AIC       BIC     logLik  deviance   Chisq Df Pr(>Chisq)
chadha.lmer2    4 374.58556 383.15810 -183.29278 366.58556                      
chadha.lmer1    5 375.75771 386.47338 -182.87885 365.75771 0.82785  1    0.36289

AICc(chadha.lmer1, chadha.lmer2)

Focus on random intercept model - refit with REML

chadha.lmer3 <- lmer(noleaves~week + (1|plant), REML=TRUE, chadha1)
summary(chadha.lmer3, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: noleaves ~ week + (1 | plant)
   Data: chadha1

REML criterion at convergence: 365.6

Scaled residuals: 
       Min         1Q     Median         3Q        Max 
-1.9226193 -0.5085957  0.0859208  0.3199276  3.9265077 

Random effects:
 Groups   Name        Variance Std.Dev.
 plant    (Intercept) 11.33340 3.366512
 Residual             16.34859 4.043338
Number of obs: 63, groups:  plant, 7

Fixed effects:
              Estimate Std. Error         df  t value   Pr(>|t|)    
(Intercept) 16.4007937  1.6886945 13.4330609  9.71211 1.9106e-07 ***
week         3.4500000  0.1972947 55.0000000 17.48653 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
week -0.584

anova(chadha.lmer3, type=3, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
      Sum Sq Mean Sq NumDF DenDF   F value     Pr(>F)    
week 4999.05 4999.05     1    55 305.77874 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get variance components

chadha.ci3 <- confint.merMod(chadha.lmer3, oldNames=FALSE)

Computing profile confidence intervals ...

chadha.vc3 <- (chadha.ci3)^2
print(chadha.vc3)

                             2.5 %       97.5 %
sd_(Intercept)|plant   2.765719795  38.24037253
sigma                 11.327647320  23.82821615
(Intercept)          170.118316638 390.40408940
week                   9.364274262  14.74475390

Random slopes output for comparison - refit using REML

chadha.lmer4 <- lmer(noleaves~week + (week||plant), REML=TRUE, chadha1)
summary(chadha.lmer4, ddf="Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: noleaves ~ week + (week || plant)
   Data: chadha1

REML criterion at convergence: 364.6

Scaled residuals: 
       Min         1Q     Median         3Q        Max 
-2.0165627 -0.5211774  0.0766027  0.3747093  3.9247059 

Random effects:
 Groups   Name        Variance   Std.Dev.
 plant    (Intercept)  6.7968482 2.607077
 plant.1  week         0.1409567 0.375442
 Residual             15.6441619 3.955270
Number of obs: 63, groups:  plant, 7

Fixed effects:
             Estimate Std. Error        df  t value   Pr(>|t|)    
(Intercept) 16.400794   1.466458 10.281494 11.18395 4.4146e-07 ***
week         3.450000   0.239551 10.281494 14.40194 3.7892e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
     (Intr)
week -0.530

anova(chadha.lmer4, type=3, ddf="Kenward-Roger")

Type III Analysis of Variance Table with Kenward-Roger's method
        Sum Sq   Mean Sq NumDF     DenDF   F value     Pr(>F)    
week 3244.8493 3244.8493     1 10.281494 207.41599 3.7892e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Check whether AR(1) covariance structure would improve the fit using nlme

chadha.lme1 <- lme(noleaves~week, random=~1|plant, method="ML", chadha1)
chadha.lme2 <- lme(noleaves~week, random=~1|plant, method="ML", correlation=corAR1(form=~1|plant), chadha1)
anova(chadha.lme1, chadha.lme2)

AICc(chadha.lme1, chadha.lme2)

no improvement with AR(1)

chadha <- lm(noleaves~week, chadha1)
summary(chadha)

Call:
lm(formula = noleaves ~ week, data = chadha1)

Residuals:
       Min         1Q     Median         3Q        Max 
-9.1007937 -3.1507937  0.2492063  2.9492063 14.5492063 

Coefficients:
              Estimate Std. Error  t value   Pr(>|t|)    
(Intercept) 16.4007937  1.4103462 11.62891 < 2.22e-16 ***
week         3.4500000  0.2506251 13.76558 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.136286 on 61 degrees of freedom
Multiple R-squared:  0.7564785, Adjusted R-squared:  0.7524863 
F-statistic: 189.4912 on 1 and 61 DF,  p-value: < 2.2204e-16